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Abstract 

Truncated Dyson-Schwinger equations represent finite subsets of the equations of 
\Q ', motion for Green's functions. Solutions to these non-linear integral equations can 

account for non-perturbative correlations. A closed set of coupled Dyson-Schwinger 
^-j- ■ equations for the propagators of gluons and ghosts in Landau gauge QCD is obtained 

by neglecting all contributions from irreducible 4-point correlations and by imple- 



menting the Slavnov-Taylor identities for the 3-point vertex functions. We solve 
this coupled set in an one-dimensional approximation which allows for an analytic 
infrared expansion necessary to obtain numerically stable results. This technique, 
which was also used in our previous solution of the gluon Dyson-Schwinger equa- 
tion in the Mandelstam approximation, is here extended to solve the coupled set of 
integral equations for the propagators of gluons and ghosts simultaneously. In partic- 
ular, the gluon propagator is shown to vanish for small spacelike momenta whereas 
the previoulsy neglected ghost propagator is found to be enhanced in the infrared. 
The running coupling of the non-perturbative subtraction scheme approaches an 



infrared stable fixed point at a critical value of the coupling, a c ~ 9.5. 
PACS Numbers: 02.30Rz 11.15.Tk 12.38.Aw 14.70.Dj 

PROGRAM SUMMARY 

Title of program: gluonghost 
Catalogue identifier: 

Program obtainable from: CPC Program Library, Queen's University of Belfast, 
N. Ireland 

Computers: Workstation DEC Alpha 500 
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Operating system under which the program has been tested: UNIX 

Programming language used: Fortran 90 

Memory required to execute with typical data: 200 kB 

No. of bits in a word: 32 

No. of processors used: 1 

Has the code been vectorized of parallelized?: No 
Peripherals used: standard output, disk 

No. of lines in distributed program, including test data, etc.: 247 

Distribution format: ASCII 

Keywords: 

Non-perturbative QCD, Dyson-Schwinger equations, gluon and ghost prop- 
agator, Landau gauge, Mandelstam approximation, non-linear integral equa- 
tions, infrared asymptotic series, constrained iterative solution. 

Nature of physical problem: 

One non-perturbative approach to non-Abelian gauge theories is to investi- 
gate their Dyson-Schwinger equations in suitable truncation schemes. For the 
pure gauge theory, i.e., for gluons and ghosts in Landau gauge QCD without 
quarks, such a scheme is derived in Ref . [1] . In numerical solutions one gener- 
ally encounters non-linear, infrared singular sets of coupled integral equations. 

Method of solution: 

The singular part of the integral equations is treated analytically and trans- 
formed into constraints extending our previous work [2] to a coupled system 
of equations. The solution in the infrared is then expanded into an asymptotic 
series which together with the known ultraviolet behavior makes a numerical 
solution tractable. 

Restrictions on the complexity of the problem: 

Solving the coupled system of Dyson-Schwinger equations relies on a modified 
angle approximation to reduce the 4-dimensional integrals to one-dimensional 
ones. 

Typical running time: One minute 



2 



References 



[1] L. von Smekal, A. Hauck and R. Alkofer, Phys. Rev. Lett. 79 (1997), 3591; 
L. von Smekal, A. Hauck and R. Alkofer, A Solution to Coupled Dyson- 



Schwinger Equations for Gluons and Ghosts in Landau Gauge, hep-ph 9707327, 
e-print, submitted to Ann. Phys., and references therein. 

[2] A. Hauck, L. von Smekal and R. Alkofer, Solving the Gluon Dyson-Schwinger 
Equation in the Mandelstam Approximation, submitted to Computer Physics 
Communications . 



LONG WRITE-UP 



1 The physical problem 



1.1 Introduction 



The infrared regime of non-Abelian gauge theories is inaccessible to perturba- 
tion theory. Confinement, being a long-distance effect, is expected to manifest 
itself in the infrared behavior of the Green's functions of the theory. Thus, 
in solving truncated sets of Dyson-Schwinger equations (DSEs) in order to 
determine the propagators self-consistently, infrared singularities have to be 
anticipated. This implies some special precautions in the numerical problem. 

In this paper we present the numerical solution of the coupled gluon and ghost 
DSEs in which the infrared behavior of the corresponding propagators is deter- 
mined analytically. In particular, asymptotic series for their infrared structure 
are calculated recursively prior to the iterative process. The DSEs being non- 
linear integral equations these series represent a systematic formulation of the 
consistency requirements in the extreme infrared on possible solutions. The 
method of simultaneously expanding the solutions to a coupled set of equa- 
tions generalizes the one used in Ref. [1] where only the gluon propagator 
was considered in an approximation in which ghost contributions are omitted 
[2-4]. 

This paper is organized as follows: In the next subsection we summarize the 
truncation scheme used in order to arrive at a closed system of equations. In 
the following subsection the reduction to one-dimensional integral equations 
is presented. For completeness we give also the most important steps of the 
renormalization procedure. In section two the numerical method is discussed 
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with special emphasis on the semi-analytic solutions in the infrared as well as 
in the ultraviolet. The numerical method based on iteration for intermediate 
momenta and matching to analytic expressions for small and large momenta 
is explained. Numerical results are presented and some implications of their 
infrared behavior are discussed. For more details in the derivation of the trun- 
cation scheme and for a further discussion of the physical implications of the 
results we refer the reader to Refs. [5,6]. 

1.2 A solvable truncation scheme 

In order to keep this paper self-contained we first summarize the trunctation 
scheme used to arrive at a closed system of equations [5]. 

For simplicity we consider the pure gauge theory and neglect all quark con- 
tributions. In addition to the elementary two-point functions, the ghost and 
gluon propagators, the Dyson-Schwinger equation for the gluon propagator 
also involves the three- and four-point vertex functions which obey their own 
Dyson-Schwinger equations. These equations involve successive higher n-point 
functions. The used truncation of the gluon equation includes to neglect all 
terms with four-gluon vertices. These are the momentum independent tad- 
pole term, an irrelevant constant which vanishes perturbatively in Landau 
gauge, and explicit two-loop contributions to the gluon DSE. The renormal- 
ized equation for the inverse gluon propagator in Euclidean momentum space 
with positive definite metric, = 5^, (color indices suppressed) is then 
given by 



D^ 1 (k) = Z 3 D tl ; 1 u (k)-g 2 N c Z 1 j* iqil D G (p) D G {q)G u {q )P ) 

+ ^N e Z i y^T% a (k,^,q)D aP (q)DMr^(-q,p,-k)^ (1) 



where p = k + q, D and T are the tree-level propagator and three-gluon 
vertex, D G is the ghost propagator and T and G are the fully dressed 3-point 
vertex functions. The DSE for the ghost propagator in Landau gauge QCD, 
without any truncations, is given by 



D G \k) = -Z 3 k 2 + g 2 N c Z, / —\ A ik» D, u (k - q) G u (k, q) D G (q) . (2) 



The coupled set of equations for the gluon and ghost propagator, eqs. (1) and 
(2), is graphically depicted in Fig. 1. The renormalized propagators for ghosts 
and gluons, D G and D, and the renormalized coupling g are defined from 
the respective bare quantities by introducing multiplicative renormalization 
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Fig. 1. Diagrammatic representation of the gluon and ghost Dyson-Schwinger equa- 
tions in the truncation scheme applied in this paper. 

constants, 

Z 3 D G :=D° G , Z 3 D^ := D% , Z g g := g . (3) 

~ 1/2 ~ 

In Landau gauge, which we adopt in the following, one has Zi = Z g Z 3 Z 3 = 1 
and Zi = ZgZ^J 2 . The SU(N C = 3) structure constants / abc of the gauge 
group (and the coupling g) are separated from the 3-point vertex functions 
by defining: 



r a ;: p (k,p,q) ^ gf abc (2n)^k+p + q)r^ p (k,p,q) . 




(4) 



The arguments of the 3-gluon vertex denote the three incoming gluon mo- 
menta according to its Lorentz indices (counter clockwise starting from the 
dot). With this definition, the tree-level vertex has the form, 

r !l P (^ P, q) = -i(k- p) p 5 flu - i(p - q)^5 vp - i(q - k)„5 w . (5) 

The arguments of the ghost-gluon vertex are the outgoing and incoming ghost 
momenta respectively, 

G^(p,q)=gf abc G,(p,q). l k_P " q (6) 



Note that the color structure of all three loop diagrams in Fig. 1 is simply 
given by f acd f bdc = —]\j c fi ab which was used in Eqs. (1) and (2) suppressing 
the trivial color structure of the propagators ~ 5 ab . 
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The ghost and gluon propagators in Landau gauge are parameterized by their 
respective renormalization functions G and Z, 



(7) 
(8) 

In order to obtain a closed set of equations for the functions G and Z the ghost- 
gluon and the 3-gluon vertex functions have to be specified. We construct 
these vertex functions from their respective Slavnov-Taylor identities (STIs) as 
entailed by the Becchi-Rouet-Stora (BRS) symmetry. In particular, in Ref. [5] 
we derive the STI of the gluon-ghost vertex from BRS invariance. Neglecting 
irreducible ghost rescattering, an assumption fully compatible with the present 
truncation scheme, this new identity together with the symmetry properties 
of the vertex constrain the full structure of the gluon-ghost vertex which 
expressed in terms of the ghost renormalization function reads [5] : 

^ / N • G ( k2 ) ■ f G ( k2 ) ,\ 

G M (p,g)=^^y+^l^py-l| . (9) 

At the same time, this gluon-ghost vertex implies a rather simple form for a 
ghost-gluon scattering kernel of tree-level structure which in turn allows for 
a straightforward resolution of the STI for the 3-gluon vertex [5]. Neglecting 
some unconstrained terms which are transverse with respect to all three gluon 
momenta the solution for the 3-gluon vertex follows from the general construc- 
tions in Refs. [7-9]. As a result, the 3-gluon vertex can again be expressed in 
terms of the gluon and ghost renormalization functions [5]: 

1W(P> 9» k ) = - A +(P 2 > <? 2 ; k 2 ) <W(P - q) P ~ A_(p 2 , q 2 ; k 2 ) 5^ i(p + q) p 
A_(p 2 , q 2 ; k 2 ) 

- 2 ^ 2 ^ (S^vpq - p u q^) i(p - q) p + cyclic permutations , (10) 

with 

A + (V 2 q 2 - k 2 ) ~ ^-1 ( G{q2) ± G{P2) \ (11) 

This establishes a closed system of equations for the renormalization func- 
tions G(k 2 ) and Z(k 2 ) of ghosts and gluons, which consists of their respective 
DSEs (1) and (2) using the vertex functions (9) and (10/11). Thereby explicit 
4-gluon vertices (in the gluon DSE (1)) as well as irreducible 4-ghost correla- 
tions (in the identity for the ghost-gluon vertex) and non-trivial contributions 
from the ghost-gluon scattering kernel (to the Slavnov-Taylor identity for the 
3-gluon vertex) were neglected. Since, at present, we do not attempt to solve 
this system in its full 4-dimensional form (but in a one-dimensional approxi- 
mation), we refer the reader to Ref. [5] for its explicit form. 



D G {k) 
D pu {k) 



G(k 2 ) 



Z(k 



k 2 
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1.3 The modified angle approximation 

In this section, illustrated at the example of the less complex ghost DSE, we 
present the approximation used to render the integral equations one-dimen- 
sional. This especially allows for a thorough discussion of their infrared and 
ultraviolet asymptotic behavior, which is a necessary prerequisite to stable 
numerical results. The leading order of the integrands in the infrared limit 
of integration momenta is hereby preserved. Furthermore, the correct short 
distance behavior of the solutions (obtained at high integration momenta) is 
also unaffected. From (2) with the vertex (9) we obtain the following equation 
for the ghost renormalization function G(k 2 ), 



where V^ip) = 5^ u — p^p v /p 2 is the transversal projector. In order to perform 
the integration over the 4-dimensional angular variables analytically, we make 
the following approximation: 

For q 2 > k 2 we assume that the functions Z and G are slowly varying with their 
arguments and that we are thus allowed to replace G{p 2 ) ~ G{k 2 ) — > G(q 2 ). 
This assumption ensures the correct leading ultraviolet behavior of the equa- 
tion according to the resummed perturbative result at one-loop level [5]. For 
all momenta being large, i.e. in the perturbative limit, this approximation 
is well justified by the slow logarithmic momentum dependence of the per- 
turbative renormalization functions for ghosts and gluons. Our solutions will 
resemble this behavior, justifying the validity of the approximation in this 
limit. Note that previously this same assumption was used for arbitrary inte- 
gration momenta q 2 in the derivation of the one-dimensional equation for the 
gluon renormalization function in Mandelstam approximation [1-4]. In this 
case, in particular for small q 2 < k 2 , the infrared enhanced solution tends to 
invalidate this assumption. 

For q 2 < k 2 we therefore proceed with an angle approximation instead, which 
preserves the limit q 2 — > of the integrands replacing the arguments of 
the functions Z and G according to G{p 2 ) = G((k — q) 2 ) — > G{k 2 ) and 
Z{p 2 ) — > Z{k 2 ). In this form, the one-dimensional approximation was used 
in a very recent investigation of the coupled system of ghost and gluon DSEs 
using only tree-level vertices [10]. However, using this approximation for ar- 
bitrary q 2 (in particular also for q 2 > k 2 as in Ref. [10]) one does not recover 
the renormalization group improved one-loop results for asymptotically large 
momenta. 




p = k-q, (12) 
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We therefore use that particular version of the two different one-dimensional 
approximations that is appropriate for the respective cases, q 2 k 2 . In this 
modified angle approximation, we obtain from (12) upon angular integration, 



1 



G(k 2 ) 



g 2 SN C l 



, (is) 

fc 2 / 



16tt 2 4 



where we introduced an 0(4)-invariant momentum cutoff A to account for 
the logarithmic ultraviolet divergence, which will have to be absorbed by the 
renormalization constant. 



It has several advantages (summarized in Ref. [5]) to use the projector 

k,,k,, 



(k) = 6^-4 



(14) 



in the gluon DSE to isolate a scalar equation for Z{k 2 ) from Eq. (1). With 
the same one-dimensional reduction as used for the ghost DSE we obtain 



i g 2 N c 

'zW)- Z3 + Zl mT 2 Y 

A 2 



k fdq 2 (7q 4 17q 2 9\ 2 2 

J v-[2¥-T¥-8) Z{q )G{q) 



+ 



fc 2 



dq*_ 7tf 



7)Z(q 2 )G(q 2 ) 



+ 



g 2 n c 

lQn 2 3 



k 2 



(15) 



In the derivation of Eq. (15), however, we omitted one contribution from the 
3-gluon loop for q 2 < k 2 , namely the following term: 



g 2 z t N c 



A Ar/ _ 2 2l 2, fZ(p 2 )G(p 2 ) Z(q 2 )G(q 2 )\ G(k 2 ) 



q 2 <k2 



N(p 2 ,q 2 ;k 2 



G(q 2 ) 



G(p 2 



k 2 p 2 q 2 



(16) 



with 



, T , , 5x 3 + 41x 2 y + 5xy 2 - 3y 3 x 2 - lOxy + 24y 2 
N(x,y;z) = — ; + 



+ 



4x(y — x) 2(y — z) 

x 3 + 9x 2 y — 9xy 2 — y 3 (2x 2 + Wxy — 3y 2 )z (x + y)z 2 



+ 



xz 



2x(x — y) 



+ 



4x(y — x) 



■ (17) 



Due to the singularity in N(p 2 , q 2 ; k 2 ) for p 2 — > q 2 which has to be cancelled 
from the terms in brackets, this contribution would generate an artificial sin- 
gularity if the angle approximation was applied. We will assess the influence 
of this term below in order to justify its omission. 
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The only difference in the 3-gluon loop as obtained here versus the Man- 
delstam approximation (see [1]) is that the gluon renormalization function 
Z is replaced by the product ZG. The system of equations (13) and (15) is 
a direct extension to the gluon DSE in the Mandelstam approximation [1]. 
Thus, the methods developed for solving the Mandelstam equation have to be 
generalized to solve the coupled Eqs. (13) and (15). 

It will furthermore become clear shortly that the leading infrared behavior of 
the solutions is unaffected by the additional manipulation to the 3-gluon loop. 
This was also confirmed in Ref. [10] where the same qualitative behavior of the 
solutions in the infrared was obtained neglecting the 3-gluon loop completely. 

With the Ansatz that for x := k 2 — > the product Z(x)G(x) — > cx K , the ghost 
DSE (13) with N c = 3 yields, 



-l 



^(^U-*)) e ~ v *- (18) 

9 /i r 



z ^\f^{n-2))^- (19) 

Furthermore, in order to obtain a positive definite function G(x) for positive 
x from a positive definite Z( — > 0, we find the necessary condition 
1/k— 1/2 > which is equivalent to 

< k < 2 . (20) 

The special case k — leads to a logarithmic singularity in Eq. (13) for x — > 0. 
In particular, assuming that ZG = c with some constant c > and x < xq for 
a sufficiently small x , we obtain G~ 1 (x) — > c (9g 2 /647r 2 ) ln(x/x ) + const and 
thus G(x) — > 0~ for rr — > 0, showing that no positive definite solution can be 
found in this case either. 

It is important to note that the ghost-loop gives infrared singular contribu- 
tions ~ x~ 2k to the gluon equation (15) while the 3-gluon loop yields terms 
proportional to x K as x — > 0, which are thus subleading contributions to the 
gluon equation in the infrared. With Eq. (18) the leading asymptotic behavior 
of Eq. (15) for x — > leads to 

2 9 9 /1 1\ 2 /3 1 1 1 V 1 2 2K 

Consistency between (19) and (21) requires that 

1 1 \ 1 9/1 1\ 

+ 77: =7 --o ■ ( 22 ) 



22-K 3 4k/ 4 2 
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Using the constraint (20) in addition, the solution is given uniquely by 

= 61-VI897 

19 v ' 

From these considerations alone we can conclude that the leading behavior of 
the gluon and ghost renormalization functions and thus their propagators in 
the infrared is entirely due to ghost contributions. The details of the treatment 
of the 3-gluon loop have no influence on the above considerations. This is 
in remarkable contrast to the Mandelstam approximation, in which the 3- 
gluon loop alone determines the infrared behavior of the gluon propagator and 
the running coupling in Landau gauge [1-4]. On the other hand, the present 
picture is confirmed by the ghost-loop only approximation to the coupled 
gluon and ghost DSEs which yields the same qualitative infrared behavior as 
investigated in Ref. [10]. The quantitative discrepancy in their numerical value 
for the exponent k ~ 0.77 can be attributed to their using of tree-level vertices 
as compared to the STI improvements used here. In contrast to the infrared, 
however, the 3-gluon loop is crucial for the correct anomalous dimensions 
which determine the leading behavior of the propagators in the ultraviolet. 

1-4 Renormalization 

In Landau gauge the renormalization constants (as introduced in Eq. (3)) obey 
the identity [11]: 

Z x = Z g Zl /2 Z 3 = 1 . (24) 

The Slavnov-Taylor identity for the ghost-gluon vertex ensures that this re- 
mains valid also in general covariant gauges as long as irreducible 4-ghost 
correlations are neglected [5]. In the following we will exploit the implication 
of Eq. (24), namely that the product g 2 Z(k 2 )G 2 (k 2 ) is renormalization group 
invariant. Near the ultraviolet fixed point this invariant is identified with the 
running coupling. Non-perturbatively, though there is no unique (scheme in- 
dependent) way of defining a running coupling, invariance under arbitrary 
renormalization group transformations (g,n) — > (g', fi') allows the identifica- 
tionQ 

g 2 Z{u,' 2 )G 2 ^' 2 ) = g' 2 = g 2 Wlp),g) . (25) 

This being one of the conditions that fix the non-perturbative subtraction 
scheme, it yields a physically sensible definition of a non-perturbative running 

1 This argument relies of course on the absence of any dimensionful parameters, 
i.e., quark masses. 
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coupling in pure Landau gauge QCD [5]. Note that this identification of the 
non-perturbative running coupling is an extension to the procedure we used in 
the Mandelstam approximation [1]. In this approximation without ghosts the 
identity Z g Z^ = 1 implies that gZ{k 2 ) is the corresponding renormalization 
group invariant product which is replaced by g 2 ZG 2 in presence of ghosts. 

The one-dimensional DSEs (13) and (15) can actually be cast in an explic- 
itly scale independent form using the following Ansatz to parameterize the 
functions G and Z motivated from their one-loop scaling behavior: 



F(x) 



1-25 

>2 




R\x) , (26) 
with x := k 2 /a and s := li 2 /a , 

where a is some currently unfixed renormalization group invariant scale pa- 
rameter and 5 = 9/44. From the definition of the running coupling (25) we 
find that g 2 (tk, g) ~ F(x) with — \ In k 2 / ji 2 . We fix the constant of propor- 
tionality for later convenience by setting (with /3 = lliV c /(487r 2 ) for Nf = 
quark flavors), 

(3 g 2 (t k ,g) = F(x) and a s (n) = £ = -j^- F(s) . (27) 

The fact that from the resulting equations besides the running coupling F(x) 
also the second function R(x) turns out to be independent on the renormaliza- 
tion scale s shows that the solutions to the renormlized DSEs for ghosts and 
gluons formally obey one-loop scaling at all scales [5]. The non-perturvative 
nature of the result thus is entirely contained in the running coupling. 

As the infrared behavior of the solutions G and Z, Eqs. (18) and (19) respec- 
tively, can be extracted without actually solving the DSEs, we find for the 
running coupling accordingly, 

g >Z(k*)G\e) = f( ti , g ) ( JLj (i - i))"' =: gj . (28) 

With Eq. (23) for k we obtain g 2 ~ 119.1 which corresponds to a critical cou- 
pling a c = g 2 /(A7i) ~ 9.48. This is in clear contrast to the running coupling 
obtained in the Mandelstam approximation [1]. The dynamical inclusion of 
ghosts changes the infrared singular coupling of the Mandelstam approxima- 
tion to an infrared finite one implying the existence of an infrared stable fixed 
point. 

With the parameterization (26) and setting k 2 = fi 2 (44> x — s) with j3 g 2 = 
F(s) we obtain equations for the renormalization constants Z s and Z 3 . For 
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the latter, this can immediately be used to eliminate Z 3 from the ghost DSE 
which then reads [5], 

^-s/f^'^-i^ 1 "^' (29) 

The gluon DSE (15) for a; = s contains the additional renormalization con- 
stant of the 3-gluon vertex Z\ which is a divergent quantity in perturbation 
theory since in Landau gauge Z\ = Z 3 /Z 3 ~ (g 2 /9o) 13S - It turns out that 
the corresponding (one-loop) renormalization scale dependence of this con- 
stant is needed in the DSE for the solution to reproduce the correct scale 
dependence asymptotically. Not so, however, a possible cutoff dependence of 
Z\ (from (?q) which cannot be removed from equation (15) consistently [5]. 
Substituting in Z\ the cutoff scale by the integration momentum y by using 
Z\ = (F (s) / F (y)) 1 ^ 36 takes care of the scale dependence of Z\ without intro- 
ducing an additional divergence. While this manipulation leads to the correct 
scaling limit for the gluon propagator [5], it might give indications towards 
possible improvements on the truncation and approximation scheme. It also 
allows to remove the gluon renormalization constant Z% from Eq. (15) and the 
same steps as for the ghost equation yield, 

11 _ jdy f7y 2 17 y 9 x\ 2S 
R 2 (x)F 1 ~ 25 (x) = J ~ \2^ ~Yx~8 + y) ^ ^ 



o 

oo 



1 l^Lvii )F 2S ( 3 fS ( x ) fdyy F s {y) l F 2S (x) 
+ X 8jy 2 ~ {V) {V) + 2li(x) J ~xl^yj 3i^T 

1 Jdy (F 2S (y) 1 a 25 \ 1 a 25 



2 J y \R 2 (y) An b 2 y 2K J An b 2 x 2K 



(30) 



o 



where k is the exponent (23). Furthermore, a := (3 g 2 = F(0) and b is defined 
through the leading infrared behavior of R(x) — > bx K for x — > 0. 



2 Numerical methods 



2.1 Asymptotic series in the infrared and behavior in the ultraviolett 



Eqs. (29) and (30) do not depend on the renormalization scale s = fi 2 /cr. This 
implies that the functions F(x) and R(x) are renormalization group invariant. 
In particular, the scaling behavior of the propagators follows trivially from the 
solution for the non-perturbative running coupling. 
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Following the method used to obtain our previous solution to the gluon DSE 
in Mandelstam approximation we are going to expand the functions F(x) and 
R(x) in the infrared in terms of asymptotic series. Due to the nature of the 
coupled set of equations a recursive calculation of the respective coefficients 
is considerably more difficult than in Mandelstam approximation [3,1]. For 
x < xq where xq is some infrared matching point, the asymptotic series to at 
least next-to-leading order is used in obtaining iterative solutions for x > xq. 
The matching point x has to be sufficiently small for the asymptotic series to 
provide the desired accuracy. On the other hand, limited by numerical stability, 
it cannot be chosen arbitrarily small either. This leads to a certain range of 
values of xq for which stable solutions are obtained with no matching point 
dependence to fixed accuracy. The additional inclusion of the next-to-next- 
to-leading order contributions in the asymptotic series has no effect other than 
increasing the allowed range for the matching point as we will show below. 

We proceed further by noting that the equation for the ghost propagator, 
Eq. (29), can be converted into a first order homogeneous linear differential 
equation for R(x) by differentiating Eq. (29) with respect to x: 

R(x ^TTW^ + F T- l -^ F ) R[x) - (31) 



The gluon equation, Eq. (30) can be rewritten as 
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R 2 (x)F 1 ~ 2S (x) 

l_x]_ 
8y 2 




17 y 9 x 7 x 2 \ , , 
- U - - - + - --)R{y)F 25 {y) 
2 x 8 y 8 y 2 I 



+ -—ba 25 y h 



7 ba 



26 



l F 2S {x) 
3 R 2 {x) 



81 
1 

" 2 



x K + A x + 



3^) fdyy F s (y) 



K 



I 



2 R(x) I x x R(y) 



,26 



,25 



dy f F 2S (y) _ J 

y \ R 2 (y) 4k b 2 y 2K J ' An b 2 x 2n ' 



+ 



(32) 



where we have used that 



x ll7 R{y)F25{y) 



_ x l f ^ ( RF *S _ ba 25 «\ + l^_ x * +Ax 

8 J y 2 V y ' 81 — k 

o y 

_ oo 

with A=^J% (RF 25 - ba 2S y K ) ■ (33) 



It follows from the leading infrared behavior, i.e., F — > a and R — > bx K for 
x — > 0, and Eq. (29) that an asymptotic infrared expansion of the l.h.s. of Eq. 
(32) has to contain powers of x K as well as integer powers of x in subsequent 
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subleading terms. This motivates the following Ansatz, 



R(x) 



bx 




,mv+'inK+l(l+2K) 



(34) 




a 



^ ' Dimn x 



,mv+'inK+l(l+2K) 



(35) 



l,m,n=0 



with Cooo = -Dooo = 1- The terms proportional to x v in these expansions 
are allowed to find the most general subleading behavior compatible with 
the consistency in the infrared. Below we will see that v ~ 2.05. Using 
2 < 3/t < 1 + 2/t < 3 one finds that different orders in this expansions do not 
mix in their successive importance at small x. Furthermore the leading infrared 
contributions are analytically evaluated and explicitly subtracted from all in- 
tegrals, assuming that the remaining contributions are integrable for x — > 0. 
For the subleading contributions to R and F suppressed by powers of x u with 
v ~ 2.05 this is justified a posteriori. 

Inserting the series (34) and (35) into Eq. (31) allows to relate the coefficients 
Ci mn to Di mn . In the solution of Eq. (31) the integration constant is set to b. 
In the order N — 1 one thus obtains: 



At higher orders in N this procedure recursively yields relations that uniquely 
determine the coefficients Ci mn in terms of the coefficients D\ mn . Further re- 
lations are obtained from Eq. (32) by expanding all ratios of R and F which 
occur with dependence on x and y, and by comparison of the respective orders, 

( x mv+(3n~2) K +l(l+2 K )j QR both gides _ Tq lead j ng Qrder N = 0{x~ 2k ), 

from Eq. (32) one obtains 



With a = p g 2 = ((9/44)(l//c — 1/2)) 1 this is nothing more than our previ- 
ously used Eq. (23) which determines k. At order N — 1 Eq. (32) yields 




(36) 
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3 1 1 1 \ a 2i 
2 2- k~ 3 + 4k J ~& 



(37) 



b 2 a l-28 
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v—1k\ 



Oix 



0{x K ) : 



0{x): 



— (Ano + 2(C 010 - 5D 010 )) = 

(X 

' + J-) -I l — ) (Quo - 5D 010 ) 



2 \2 + z/ - re 2 — re/ 3 v - 2k 

— (D 001 + 2(C 00 i - SD 001 )) = -b 3 f( K ) + 



1 



X - ) (Cooi - SD Q01 ) 



2 V2 + 2re 2-kJ 3 k 

— (D 100 + 2(C 100 - 8D 100 )) = -\tA + 



1 



2 \3 + k 2-k 



2 ) (Cioo — 3D 100 ) , 



(38) 



(39) 



(40) 



with 



/(«) == 
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7 

+ - + 



2(3 + k) 2(2 + k) 8(1 + k) re 8(1 - re) ' 



(41) 



Eqs. (36) and (38)-(41) determine the coefficients Z) and C to lowest non- 
trivial order. They decouple into three times two equations for each pair 
(C lmn ,D lmn ). For (l,m,n) = (1,0,0) one obtains 

C 100 ^ 0.055546 2 A and D wo ~ -0.69926 2 A , (42) 

where the constant A is defined in Eq. (33) The set of equations for (/, m, n) = 
(0, 1,0), Eqs. (38)-(41), is homogeneous. The determinant of its 2-dimensional 
coefficient matrix is zero for 



-6 - re - 3re 2 ± ^(3 + 2re) (104 + 92re) re 2 + (6 + re + 3re 2 ) 2 

2 (3 + 2re) ' ^ 3 ^ 

There exists one positive root for the plus sign which determines the positive 
exponent v. With re = (61 — v / 1897)/19 one arrives at 

v ~ 2.051 and C ow = -0.0124D io • (44) 

For (/, m, n) = (0, 0, 1) the scale of the coefficients is set by the inhomogeneity 
in Eq. (39), b 3 f(n), yielding 

Cooi ^ 1-969 b 3 , and D 001 ~ -26.52 b 3 . (45) 

Based on these next -to-leading order results, higher orders, though increas- 
ingly tedious, can be obtained recursively by analogous sets of equations. The 
general pattern is such that the lower order fixes the scales for higher or- 
der coefficients. This allows to define scale independent coefficients C and 
D by extracting their respective scales according to the exponent Ti mn = 
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mv + SnK + 1(1 + 2k) of x for a given set (I, m, n), 

C lmn =:C lmn b 3n+2l t m A l , and D lmn =: D lmn b 3n+21 t m A 1 . (46) 

The scale of the powers of x v is set by D 010 and we have defined for convenience 

* := -Ano (47) 

i.e., Coio — 0.0124 and D ow = —1. We summarize the values of the coefficients 
C and D for JV = 2 in table 2.1. 



(Z,m,n) 


(2,0,0) 


(1,1,0) 


(1,0,1) 


(0,2,0) (0,1,1) 


(0,0,2) 


C 


-0.1042 


-0.3034 


-7.933 


-0.2160 -11.55 


-151.0 


5 


0.5246 


1.590 


40.10 


1.226 60.98 


766.8 



Table 1 

Coefficients of the asymptotic expansion for N = 2. 

The two parameters 6 and t are related to the overall momentum scale. After 
fixing the scale one is still left with one independent parameter. This leads to a 
scale invariance which can be described as follows: A change in the momentum 
scale a (introduced by x = k 2 /a) according to a — > a' = a/A or, equivalently, 
x — > x' = \x can be compensated by 

b -> b' = b/X K , and t -> t' = t/X" . (48) 

We choose the scale without loss of generality such that the positive number 
6=1. The parameter t can in principle be any real number including zero. We 
can find numerically stable iterative solutions for not too large absolute values 
of t (see below) . Furthermore, it can be verified numerically, that a solution for 
a value of 6 ^ 1 for fixed t is identical to a solution for 6 = 1 and t' — t b v l K ', 
if x is substituted by x' = x/b 1 ^. This is the numerical manifestation of 
the scale invariance mentioned above (for A = 1/6 1 / K ). Note that under scale 
transformations (48) the constant A appearing in Eq. (33) trivially transforms 
according to its dimension, A — > A' = A/X, without any adjustments from 
the way it is calculated: 



A' = A/X= lim - 



8 



[%R(y')F 2S (y')-b'a 2S ^^\ . (49) 

J y l — k j 



V 



In Fig. 2 the numerical solutions for F(x) and R(x) for 6 = 1 and t — at 
small a; together with their respective asymptotic forms to order N — 1 and 
N = 2 are displayed. The contributions of the order iV = 2 in the asymptotic 
expansion become comparable in size to the lower order at about x ~ 0.2. As 
the error in the asymptotic series is of the order of the first terms neglected, 
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0(N=1) 
0(N=2) 



0.1 



1 




0.01 







0.1 



0.2 



0.3 



0.4 



0.5 



0.6 



0.7 



0.8 



Fig. 2. The numerical solutions of F(x) and R(x) for t = and 6=1 together with 
their asymptotic expansions to order N = 1 as well as N = 2 at small x. 

this supplies an estimate for the range of x in which the asymptotic expansion 
can yield reliable results. In the particular calculation described below we used 
a value of about Xo = 0.01 for the matching point relating the result of the 
iterative process to the asymptotic expansion. This is obviously well below the 
estimated range of the validity of the asymptotic expansion. 

As already stated, for intermediate momenta the integral is done numerically. 
In the ultraviolet limit, i.e. for x — > oo, we have F(x) — > 1/ \nx and R(x) — > 1 
which allows us to alleviate the cut-off dependence in the numerical determi- 
nation of A. Note that this is the only integral left with the upper boundary 
being infinity. In Eq. (33) the corresponding integral is calculated using an 
ultraviolet matching point X\ and 



for sufficiently large x±, where T(a,x) is the incomplete gamma function. 

Similarly to the case of Mandelstam's approximation [1] we calculated all 
integrals using a Simpson integration routine of fourth order. In order to reduce 
the numerical errors which can otherwise destroy the convergence we had to 
use meshpoints equidistant on a logarithmic scale, i.e., 



Convergence properties are furthemore significantly improved by weighting 
the iteration: Instead of a full update of the functions with every step we 




(50) 




(51) 
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introduced an exponentially distributed weight between the two functions, 

rj = I e -(A.+A«)/2 ; (52) 

where is defined by 

A F := max{|F m (j)/F(j) - 1|}, , (53) 

and with an analogous definition of Ar. Here, F i+i is the preliminary result 
of the ith iteration step calcualted using the input Fj as obtained from the 
previous iteration. From this, the input for the subsequent iteration is choosen 
not to be the full new F i+1 but rather 

F i+1 = (1 - rj)F i+1 + v Fi , and R i+1 = (1 - 77)^+1 + rjRi (54) 

analogously. This increases the stability of the algorithm by suppressing pos- 
sible oszillations in the iteration. 



2.2 Numerical results 



Most of the numerical results reported here were obtained with the order 
N — 1 in the asymptotic expansion. We checked explicitely for all cases that 
no dependence on the matching point exists for 0.01 < x < 0.1 We have 
calculated F and R for several values in the range — 5 < t < 16. At lower 
negative values the procedure became numerically unstable due to a develop- 
ing (tachyonic) pole in F(x). The fact that the integral equations for R and 
F possess a one-parameter family of solutions characterized by t is in fact 
the reason for the necessity of the infrared expansion up to next to leading 
order. No stable solution can be found numerically without fixing the leading 
x-dependence of F(x) at small x by choosing a value for the parameter t. 
This is a boundary condition to be imposed on the solutions from physical 
arguments. 

In Fig. 3 the numerical results are displayed for different values of the pa- 
rameter t (all with 6 = 1). Perturbatively, we expect R(x) to approach a 
constant value and F(x) — > l/ln(Xx) for x — > 00. The reason we introduced 
here the constant A in the one-loop running coupling F is that we fixed the 
momentum scale in our calculations by arbitrarily setting 6 = 1. The relation 
between the scale of perturbative QCD Aqcd an d a cannot be determined 
this way. Therefore, we set Aq CD = <j/\ for some scale parameter A. Fixing 
the scale in our calculations from the phenomenological value of 0:5 at the 
mass of the Z-boson, one obtains a ~ (350MeV) 2 for the t — solution for 
F. A detailed discussion of the anomalous dimensions of gluons and ghosts 
allows an estimate of A to be in the range 1.5 ~ 2 which corresponds to a 
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0.1 





t = -4, -2, -1,0, 1,2, 4, 8, 16 


F(x) 


\ 


t = -4 




t = 16 






1 R(x) 


t = 16 




t = -4 



0.1 1 10 100 10 4 10 6 10 



x = k2, 

Fig. 3. The numerical solutions of F(x) and R(x) with 6=1 for different values of 
the parameter t = {—4, —2, —1, 0, 1, 2, 4, 8, 16} (solid lines represent t = solutions). 

Aqcd in 250 ~ 300MeV [5]. The solutions for t ^ display a qualitatively 
similar behavior at high momenta with slightly different values. The solutions 
for positive values of t seem to have more residual momentum dependence in 
R at high momenta than those for t < 0. For negative values of t the run- 
ning coupling, « s (/i) = F(s)/(47r/3 ), has a maximum, a max > a C) at a finite 
value of the renormalization scale s = fi 2 /a. This is because the dominant 
subleading term of the running coupling in the infrared is determined by t, 

F(x) -> a(l - tx v + D 001 x 3k + D 100 x 1+2k ) , x -> . (55) 

With v ~ 2.05 < 3k < 1 + 2k and £> oi < 0, it is clear that for t < the 
running coupling increases for smallest scales close to /i — before higher 
order terms dominate. There necessarily has to be a maximum a max > « c at 
some finite scale \i for any solution with t < 0. 

For t > 0, a c = a(/i = 0) is the only maximum of the running coupling 
for all real values of the renormalization scale, and a c is thus a true infrared 
stable fixed point. Comparing the behavior of the resulting gluon and ghost 
renormalization functions in the ultraviolet we observe that, for the t > so- 
lutions, the case t = yields the best resemblance of their one-loop anomalous 
dimensions. We therefore interpret the case t = as the most physical one and 
conceive the existence of solutions for t ^ as a mathematical peculiarity. 

In reducing the DSE for the gluon propagator to a one-dimensional equation 
we had to dismiss the contribution (16) in order to avoid an artificial sin- 
gular contribution. To asses whether this is justified we calculate the contri- 
bution from (16) without the one-dimensional approximation from the gluon 
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_o I I I I I 

10 10 10 10 10 

Fig. 4. The dismissed contribution (16) compared to the inverse gluon function. 

and ghost renormalization functions as obtained from the iterative scheme, 
i.e., from the one-dimensional integral equations. In Fig. 4 the inverse gluon 
function is displayed as a measure of the terms retained on the r.h.s. of the 
one-dimensional equation (15). This is to be compared to the neglected contri- 
bution (16) as calculated from the selfconsistent results. One clearly observes 
that the dismissed contribution remains small at all momenta and becomes 
negligible for small and large momenta quickly. Although, in principle, even 
small contributions might become important in non-linear self-consistency 
problems this is rather convincing support for the omission of the terms in 
(16) to which the one-dimensional approximation cannot be applied. 



3 Conclusions 



We have solved a set of two coupled non-linear integral equations. These 
solutions involve functions which are highly singular in the infrared. The cor- 
responding infrared behavior has to be treated analytically by converting the 
integral equations into recursion relations for the coefficients of asymptotic 
expansions. The final numerical solution is obtained by matching the asymp- 
totic expansions to an iteration process used for momenta above the matching 
point. 

The numerical algorithm described here was used in the solution of the cou- 
pled gluon-ghost Dyson-Schwinger equations reported in Refs. [5,6] for the 
first time. The resulting gluon and ghost propagators displayed a new type of 
infrared behavior involving irrational exponents of the momenta. This generic 
type of DSE solutions for gluon and ghost propagators of the same qualitative 
form has been verified recently using a different truncation scheme and differ- 
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ent numerical methods [10]. We therefore believe that the algorithm presented 
here will prove useful in further studies of Dyson-Schwinger equations also. 



4 Description of the program 

4-1 The main program 

After defining the variables and setting the parameters k, u, a etc. and the 
matching points x and x\ the functions F(x) and R(x) are initialized to 

Fix) = — — ^ and R(x) = 1 + (x - l)e~ x . (56) 

ln(l.l + x) 

The iteration process consists of several parts. The first is the evaluation of 
the constant A (see (33) and (50)) using the functions determined in the pre- 
vious iteration step. Hereby Eq. (50) is used for large momenta. Next, the 
contributions to the gluon and ghost equations due to the infrared region is 
calculated using the expansion in the asymptotic series. In the intermediate 
momentum range the integrals above the infrared matching point are com- 
puted numerically with the help of an Simpson routine of fourth order using 
the mesh defined by Eq. (51). Application of Eqs. (52) to (54) completes the 
iteration step. 

Convergence is tested by comparing the input and output functions of an 
iteration step pointwise. If the maximum relative deviation is less than EPS 
it is assumed that convergence is achieved, and the result is written to the file 
gluonghost.out in three-column form: x, F(x), R(x). 

4-2 Subroutines and functions 
Function Simpson 

Returns the integral of a function which is given at equally spaced abscissas. As 
far as the number of abscissas is sufficient this function uses a closed Simpson 
rule of order l/N 4 [13]. 

Function Gamma 

This routine returns the incomplete gamma function T(a, x) using a continued 
fraction as described in [13]. 
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5 Testing the program 



Naturally, trivial tests establishing the independence of the number of mesh- 
points, the infrared matching point xq, the ultraviolet cut-off x\ and the order 
of the asymptotic series in the infrared have been performed. We could also 
verify that the results are independent of the initializing functions chosen at 
the beginning. 
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TEST RUN 

standard output 

Number of meshpoints : 500 
Infrared matching point : 0.01 
Ultraviolet cut-off : 1.0E+08 
parameter t = 0.0 
eps = 1.0E-07 



Convergence achieved after 126 iterations! 
max. deviation between Fin and Fout : 9.47201E-08 
max. deviation between Rin and Rout: 3.42583E-08 
Output written to gluonghost . out 



0.1000000000E-01 
0.1047128548E-01 
0.1096478196E-01 

0.9120108394E+08 
0.9549925860E+08 
0.1000000000E+09 



gluonghost. out 

0.8298109605E+01 
0.8297891961E+01 
0.8298096464E+01 

0.6234215822E-01 
0.6218344927E-01 
0.6202552418E-01 



0.1457629245E-01 
0.1520601418E-01 
0.1586253844E-01 

. 9295472468E+00 
0.9296193060E+00 
0.9296910980E+00 
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